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Summary 


First- and second-order numerical proccduics for calculating two- 
dimensional transonic flows that treat shock wav*.s as disconti- 
nuities are dine. ’scd. This short commun’ cati on illustrates their 
application to a simple but non-trivial problem for which there are 
^limited theoretical results. 

Introduction 

Some of the numerical difficulties that arise in inviscid transonic 
flow computations occur because the solution that is being sought 
is discontinuous and the numerical procedures often employed ap- 
proximate the discontinuities by steep gradients; thus relatively 
small mesh sizes are required if dissipative and dispersive trunca- 
tion errors are to be acceptadsle. Further, the difference 
equations used, as Murman has pointed out [1], must be chosen with 
•care if entities that are conserved by the basic equations across 
discontinuities are to be conserved by the finite difference scheme. 

We believe there are important advantages to computing such flows 
with numerical procedures that treat discontinuities as such, pro- 
vided this can be done without major programming complexity. 

Moretti [2] has been the most constant advocate of doing precisely 
this, and he has computed a niunber of steady and unsteady flows 
using procedures . are essentially tailored to hyperbolic 
problems. Salas [3j has recently computed supersonic duct flows 
using somewhat analogous techniques, and Yu and Seebass [4] ex- 
amined a transonic flow problem where the discontinuity could be 
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traced from a boundary where the flow was hyperbolic. 




The use of similar procedures for mixed flows, with a minimum of 
additional computational complexity has been pursued by Hafez aiid 


Cheng (5], and by \is. Some of our results arc the subject of this 
brief communication. 


Transonic Fl ow Past a Wedge 

Wo have applied various numerical procedures to calculate flows 
past simple two-dimensional lifting airfoils; we illustrate one of 
them here for a very simple, but non-trivial, problem that models 
slightly subsonic flow past a wedge. The mathematical problem as- 
sociated with the small perturbation ai^proximation has been solved 
analytically for a sonic free stream, and we know the first deriv- 
ative of the drag with Mach number evaluated at sonic conditor.s 
(Guderley and Yoshihara (6]; Liepmann and Bryson t7]). With the 
use of modern computers to evaluate the hypergeometric function, it 
would seem that precise results could be obtained for a range of 
Mach numbers. 

The correct hodograph formulation of this problem was tackled 
numerically by Yoshihara [8] nearly two decades ago. v.Tiile these 
results may not be definitive enough to validate the computation of 
the supersonic domain, they are accurate in their predictions of 
the pressure coefficient on the nose of the wedge, and consequently, 
the drag. 

The mathematical formulation of this problem is simple and well- 
known. We need to solve 

(K + “ 0 subject to 

( 1 ) 

<t»^(C,0) = 0 for 5 ^ {0,1], \.K,0) =1 for ^ e [0,1]. 

The solution must behave as 

4)^^ (^,n)dCdn 

2 2 

for 5 +n ® . This expression is useful for supplying boundary 
data on a finite domain used in nvimerical computation. For small 
values of <, i.e., for near sonic flow, a less severe asymptotic 
representation must be used. As < -* 0 we can use the results of 
Reference 5. Here 4> is proportional to the perturbation velocity 


<|) 'V 


2Ti/r(^^ - <n^) 




potential, K is the usual tx^insonic similarity parameter 
2 2/3 

(M^ -1) / [(Y+1) 6j , 6 the wedge half angle and 4 ®^d n the 
usual non-dimensional coordinates. The jumo relations that follow 
directly from (1) and the irrotationality condition, are 

[K 4^ -i^)^ - <VV^' 

1 ( 2 ) 
dn/dC “ + [< + J 2 

Here denotes values downstream of the shock wave. As Oswatitsch 
[91 noted many years ago, the wave drag of a body is directly re- 
lated to the entropy produced by these waves. For this problem, 
the normalized drag coefficient may be determined, to lowest order, 
by the alternate expressions 


Cp » 4>(1,0) - <P(0,0) 


(a^6U) 


((fi^ - tf)^ ) ^(30 


where is the critical speed, U the free stream velocity, and da 
an element of shock surface. For ic •+■ 0, 1.75 + 2<. 

Numerical Procedures 


For these calculations we have used a first-order accura\,e conser- 
vative scheme to calculate the solution to Equation (1) . In 
conjunction with this scheme we have also used shock fitting. 

Because we anticipated treating the shock as a discontinuity and 
satisfying the jump relations across it, we were not initially con- 
cerned with using conservative differencing which has been shown to 
be essential across shock waves [1]. Our results indicated, however, 
that our shock-fitting procedure will not converge to the correct 
result if it is introduced into a converged non-conservative 
calculation. 


Equation (1) is of mixed type, and it is either hyperbolic or 
elliptic depending on whether tc + <})^ is positive or negative. 

Type dependent difference schemes, first introduced by Murman and 
Cole [10}, have proved effective in solving such equations. To de- 
termine the type of the equation K 4 is computed both by the 
central and by backward difference approximations. These results 


are then used to select prop<.>r difference approxinations for the 
derivatives. The only difference between the present scheme and 
the fully conservative scheme (1] is in the choice of the sonic 
point approximation. The fully conservative scheme failed because 
of the poor approximation of the differential e.juation at the 
sonic point caused by insisting upon flux conservation. Along the 
first characteristic of the expansion fan the first term in the 
differential equation (1) must precisely balance the second term 
as ic + 0. Consequently, we have found it necessary to adopt 

a non-conservative approximation of the sonic point. It is obvious, 
of course, that at the shock point the flux conservative fona must 
be used if the shock is not treated as a discontinuity. There the 
differential equation is of little consequence, but flux conserva- 
tion is. 

Shock Fitting 

When shock waves are embedded in sharp gradients, fairly refined 
mesh spacing may be required if we want accurate results from a 
flux conservative difference approximation. An alternative, as we 
have mentioned, is to treat the sliock as a discontinuity in the 
later stages of such calculations. For some problems this may prove 
to be essential; for others, it may be desirable, and for many 
simply superfluous. It may prove particularly beneficial for shock 
waves that originate with an elliptic behavior downstream and then 
undergo a transition to hyperbolic behavior there. Shock fitting 
also provides an easy mechanism for calculating the lowest order 
wave drag by integrating the entropy rise across the shock. For 
flows in which the flow behind the shock is elliptic, flux con- 
servative calculations provide a reasonable definition of the shock 
and also can be used to calculate the entropy rise. For these 
calculations we have not introduced our shock-fitting procedures 
until we have obtained a reasonadaly well-converged first-order 
result. The initial shock position is determined by interpolating 
the sonic position in the supersonic to subsonic transitioiiS, and 
the shock point (s) are then treated as regular computational 
points (see sketch) . The flow properties ahead of the shock can be 
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determined either by extra- 
polation from the upstream 
conditions at Q and R, or by 
using the characteristic re- 
lations along C ±. For those 
calculations the simple extra- 
polation method was used. 
Behind the shock 6 is ex- 

n 

trapolated by the usual dif- 


ference method using the old 

A 

value at P anc.' is obtained 
from the jump relation (2) . 


Should the shock slope extra- 



polation provide a shock intersection at s’, as indicated by the 


dotted line, then the value of ^ is fixed at s ' by using the values 

A ^ ^ 

of 4> “ 4*» and at s. At each stage of the calculation we 

correct the position of the shock by the simple procedure 



(3) 


Various values of k have been tried; values near one seem to be 
the most satisfactory. The rationale for (3) is simple: if the 

flow ahead of the shock has increased in speed, i.e., 

> 0, (3) allows the shock to be swept downstream and 
flow properties are re-calculated using the new shock position. 

Such corrections are repeated with each iteration until the changes 
in 4> are judged sufficiently small. We examined a simple one- 
dimensional problem (discussed in Bauer et al. tllj) using the pro- 
cedure, described above. We found that the shock approaches the 
exact location independent of the initial guess for its location. 

The computation time for this simple model was comparable to that of 
an equivalent shock point operator scheme, and the accuracy slightly 
better. 


Hafez and Cheng 18] have pursued somewl.at similar calculations. 

They calculated the transonic flow over a parabolic-arc airfoil 
using a shcck-fitting scheme. In their scheme the shock is located 



by a numerical procedure that detects a jump in and the dif- 
ference approximation at the point following the shock is construct- 
ed using the jump relation. Such a scheme may be inherently less 
accurate than oursj it is also loss complex and may be less sensi- 
tive to the initial data used. The '•esults they report are in good 
agreement with the fully conservative results of Murman tl] and with 
our results. If we start with a converged flux conservative calcu- 
lation and introduce shock fitting we find minor changes in the 
shock shape and the local shock slope. For flow past curved pro- 
files we find a behe.vior that closely approximates the Oswatitsch- 
Zierep (12) singularity expected behind the snock. 


Results and Discussion 


The supersonic portion of the flow field calculated using a first- 
order shock-fitting scheme is depicted in Figure 1 for < = - 0.5 
and a mesh spacing that corresponds, to twenty points on the wedge 
nose. The total computational region was “2 < ^ < 4 and 0 < D < 6 
A weak shock seems to originate near the wedge corner but it has 
been smeared out by truncation errors. There is further evidence of 

such a shock in the second-order calculations which we have carried 

* 

out. Figure 2 depicts the pressure coefficient = -2 (k + 

for three values of D. Figure 3 depicts the shock image in the 

* _ * 

(q I 0) -plane were q = k + and 0 is the normalized local flow 
deflection angle 9/6. Near the tip of the shock, the shock is of 
the supersonic-supersonic type. 


The normalized drag coefficient for K = -0.5 has been calculated 
using various schemes. The shock-fitting algorithm gives a normali- 
zed drag coefficient of 0.656 when evaluated by pressure integral, 
and 0.561 when evaluated by entropy production. The values of 
for non-conservative first- and second-order approximations are 
0.677 and 0.723 respectively. Until we are 2d3le to carry out 

calculations with more refined mesh spacings and at smaller values 

% 

of < it nuikes little sense to compare them to the result 'V' 1.75 
+ 2K for IC 0. Nevertheless, it is probably reasonable to conclude 
that this theoretical value underestimates the drag (see (8J), and 
hence so do the calculations we have reported here. 


As noted earlier, we have applied the algorithm outlined here to 
flows past airfoils. The small perturbation velocity distribution 
for a flux conservative calculation is compared in Figure 4 to that 
obtained using shock fitting for the same computational grid. 

Aside from the better definition of the singula.* behavior of the 
pressure gradient behind the shock and minor changes in shock shape 
(not shown), the results are comparable. 

The research was supported by the NASA on Grant NGR 33-010-057. 
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Fig. 4. Velocity distribution 
along the chord of a 6% thick 
parabolic arc for = 0.909 
with and without shock fitting. 




Fig. 3. Shock polar showing 
supersonic conditions behind 
tip of the shock. 
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